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Abstract 



We analyze a non-symmetric coupling of interior penalty discontinuous Galerkin 
and boundary element methods in two and three dimensions. Main results are 
' discrete coercivity of the method, and thus unique solvability, and quasi-optimal 

convergence. The proof of coercivity is based on a localized variant of the varia- 
tional technique from [F.-J. Sayas, The validity of Johnson-Nedelec's BEM-FEM 
coupling on polygonal interfaces, SIAM J. Numer. Anal, 47(5) :345 1-3463, 2009]. 
This localization gives rise to terms which are carefully analyzed in fractional order 
Sobolev spaces, and by using scaling arguments for rigid transformations. Numeri- 
■^j- | cal evidence of the proven convergence properties has been published previously. 

m ■ AMS Subject classification. 65N30, 65N38, 65N12, 65N15 

(N , 

1 Introduction 

In a recent article, Of, Rodin, Steinbach and Taus [19] propose three discretization meth- 
ods that combine Interior Penalty Discontinuous Galerkin methods (IPDG) with Bound- 
ary Element Methods (BEM). One of the methods falls into the category of non-symmetric 
coupling of Finite and Boundary Elements, while the other two belong to the general sym- 
metric coupling philosophy. Only one of the symmetric methods is analyzed but numerical 
evidence of the good properties of the non-symmetric coupling is given. In this paper we 
prove that the non-symmetric method in [19J converges, this being, to the best of our 
knowledge, the first successful analysis of a non-symmetric coupling of DG and BEM. 

Let us first briefly revise the milestones of the literature of BEM-FEM coupled schemes. 
The mathematical literature on BEM-FEM coupling can be traced back to the seminal 
work of Brezzi, Johnson and Nedelec [5j [16], the article [16] being an early and very 
relevant contribution to the subject. (It has to be noted, though, that the engineering lit- 
erature had previously visited these ideas and produced interesting results.) The method 
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of [16] uses one integral equation -derived from Green's Third Identity- to construct a 
non-local boundary condition in order to cut off the computational domain for an exterior 
diffusion problem. Because of the way the analysis was approached, using compactness 
arguments, the cut-off boundary where the non-local integral condition is imposed had 
to be taken to be smooth in order for the analytical arguments to be meaningful. While 
practitioners never found clear reasons to dismiss the use of simpler polygonal interfaces, 
theory did not move very much from this initial stagnation for more than two decades. 
Instead, symmetric coupling methods overcame the theoretical difficulty by using a sec- 
ond integral equation in the set, either mixing the integral operators with the interior 
formulation (two-field method) or by using an additional boundary unknown (three-field 
method). This restored symmetry to the coupled formulation and, with it, coercivity. The 
original work of Costabel and Stephan [9j [10] , and Han [15] started this fruitful trend that 
for many years enjoyed the prestige of full theoretical justification. In particular, most 
coupling schemes for the Boundary Element method with other domain methods (Mixed 
Finite Elements, Non- conforming FEM or Discontinuous Galerkin approximations) were 
directly based on the ideas of symmetric coupling. The monograph [13] collects much of 
what was known in the mid-nineties on the mathematics of BEM-FEM coupling. 

Non-symmetric coupling was revisited recently with the very simple result that varia- 
tional techniques were enough to prove stability of Galerkin methods for non-symmetric 
BEM-FEM coupling. The first result in this direction [20] was further simplified in [21] 
and [IT], and right now, it is clear that the proof itself does not contain any theoretical 
ingredient that does not belong to the analytical toolbox that is employed for analysis 
of elliptic boundary integral equations [17]. Similar ideas are treated in [18] for non- 
symmetric coupling of Mixed-FEM with BEM. As already mentioned, this work is the 
first contribution to non-symmetric coupling of DG and BEM. 

Mathematical theory for the coupling of Discontinuous Galerkin methods (DG) and 
BEM is less than one decade old. The coupling of Locally DG (LDG) methods with BEM 
was proposed in [T3], and extended in [6], based on a symmetric three-field formulation 
with an additional mortar variable. Although the LDG method has the aspect of a mixed 
method (it approximates both the potential and its gradient), it can be described using 
a so-called primal form (that uses only the potential). Barring technical difficulties, this 
meant that the symmetric coupling led to an analysis based on an energy (coercivity) 
estimate in the proper discrete norm. The primal formulation in [HJ [6] is non-consistent 
and a Strang-type analysis is needed. A later paper [12] eliminated the need of the mor- 
tar variable and many of the theoretical difficulties by demanding that the discontinuous 
piecewise polynomial functions that approximate the potential in the LDG method be- 
come continuous at the coupling interface. In practice this can be enforced using Lagrange 
multipliers (see some further explanations in [IS])- The paper [12] also showed how to 
generalize to some methods of the IPDG class. 

Insisting in the symmetric approach, the paper [S] proposes a systematic approach to 
couple BEM with DG-FEM and shows that most known methods fit into a double general 
framework. The symmetric coupling of Hybridizable DG (HDG) with BEM is proposed 
in [8], although analysis was postponed to the more recent [7J. Note that, in the very 
different context of the transient wave equation, and including only a stability analysis 
based on energy arguments, Abboud, Joly, Rodriguez and Terrase [I] present the first 
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time-domain coupling of DG methods with BEM. 

Coming back to the idea of non-symmetric methods, in this paper we analyze a non- 
symmetric coupling of IPDG with BEM proposed in [19J. The model problem will be a 
transmission problem in free space (in two and three dimensions). Without giving full 
details at this moment, let us explain what the difficulties are. The discrete scheme is the 
search for (u^, Xh) 6 14 x A/, such that: 

ana{uh, v h ) - (X h , v h )v = {f,v h )n_ + (Pi,v h } T Vv h eV h , ^ _ 

_ (fj-h, \uh - Ku h ) v + (/i h ,VX h ) T = (fj, h , \fh - K.(3 )t V/i h e A h . 

Here Vh is a space of discontinuous piecewise polynomial functions on a triangulation 
of the domain f2_ and is a space of piecewise polynomials of the same degree on 
the triangulation of the interface T that is inherited from the one of The bilinear 
form avG(uh,Vh) approximates the Dirichlet form (Vw, Vf)^_ and contains a first group 
of non-conforming terms, including jumps of the solution and the test function and a 
penalization term on the element faces included to stabilize the method. The class of 
methods that we treat in this general DG scheme include the original Interior Penalty 
method of Arnold j2] and some symmetric and non-symmetric variants, all of them fitting 
in the successful unified framework developed by Arnold, Brezzi, Cockburn and Marini in 
[3J. The symbols /C and V correspond to two boundary integral operators and create the 
non-local integral boundary condition for the problem. The off-diagonal terms in (II. II) . 
involving boundary and interior quantities, are actually non- conforming approximations 
of H~ l l 2 (T) x i/ 1 / 2 (r) duality products in the variational formulation of the problem. 

As already intuited in a remark in [19], the variational technique of [20] (essentially 
integration by parts after recognizing that the exterior potential has to be transmitted to 
the interior domain) plays a key role in this analysis. However, integration by parts has 
now to be applied element-by-element and a whole new array of terms have to be bounded 
below and hidden by carefully tuned weighted Young inequalities in the essential coercivity 
estimate (given in Theorem 12. ip . The analysis requires the handling of fractional order 
Sobolev spaces on the interior faces (scalability of these norms -or lack thereof- will be a 
fact to keep in mind) and some scaling arguments that have to be dealt with rigidly, due 
to the fact that harmonicity of the exterior potential that is transmitted to the interior 
domain is needed in key steps of the process. A new source of theoretical complications 
stems from the fact that the discrete norm for which the coercivity result holds does not 
contain any term involving the value of Uh on the interface V . This fact will make the 
bilinear form not bounded in the coercivity norm and a stronger norm has to be produced 
in upper bounds for the global bilinear form. A final detail that is not entirely obvious 
arises from the fact that we are dealing with a transmission problem in free space and 
that constants in the interior domain play a certain separate role from the analytical point 
of view. While this difficulty could be easily removed by considering a simpler problem 
(an exterior Dirichlet problem), we find it worthwhile to work out all the details for this 
case. The entire analysis (especially discrete coercivity) requires the careful handling of 
inequalities related to fractional order Sobolev norms in some reference configurations. 
We have striven to make all details as transparent as possible to the reader, in the hope 
that they will help for possible future generalizations. 

The paper is structured as follows. In Section [2] we present the model problem, its 
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variational formulation as a boundary-field problem and its discretization with an IPDG- 
BEM scheme. We also state the two main results of this paper: discrete coercivity 
(Theorem I2.ip and optimal convergence (Theorem 12.31) . Section [3] is devoted to proving 
Theorem I2.1[ building up from technical estimates to a full detailed proof of the discrete 
coercivity of the method. Section H] works in a very similar way to provide a proof of 
Theorem 12.31 We remark that numerical evidence of the performance of this method has 
already been presented in 1 1!) . 

Prerequisites. Some basic knowledge of the basic Sobolev spaces H m (0), their norms 
|| • \\ m> o and seminorms | ■ \ m p will be assumed throughout. A simple subscripted norm 
| ■ || B will always refer to the L 2 (B) norm. The fractional order Sobolev spaces if ±1/2 (r) 
will be used from the beginning of the paper in order to introduce the formulation and to 
present the results. However, in Section [3] we will be very precise in positive and negative 
order Sobolev spaces both on a domain, its boundary or part of it. A very detailed 
reference for these results is the monograph of McLean [17] that also includes proofs of all 
the mapping properties of the potentials and integral operators that will be loosely used 
in this work. 



2 A non-symmetric coupling of IPDG and BEM 

2.1 Model problem and two-field non-symmetric formulation 

Let Q- be a bounded polygonal domain in the plane or a polyhedral domain with Lipschitz 
boundary in the space. Let T := dfl_ and fl + := For the sake of simplicity, we will 

assume that fl + is connected. The symbol 7 will be used to denote the trace operator, 
while d v will be used for the normal derivative. The model problem is a transmission 
problem of the form 

—Au = f in 

ju = 7« + + p on T, 

d u u = d u u + + Pi on T, (2.1) 
—Au + = in Q + , 

u = 0(1/ r) as r 00. 

We assume that (3 £ ff 1 / 2 (r), f} x £ L 2 (T) and / £ L 2 (VL_). A necessary and sufficient 
condition for existence of solution in the two-dimensional case is 

/ / + / Pi = 0. (2.2) 

In the two-dimensional case, it follows from well-known results on potential theory that 
d v u + £ H 1/2 (T) := {A £ H- X '\Y) : (A, l) r = 0}, (2.3) 
where the angled brackets are used to denote the ff _1 / 2 (r) x ff 1 / 2 (r) duality product. 
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For representation of the exterior solution we need to introduce the layer potentials 
on T and some associated boundary integral operators. Let 

{-^log|x-y|, if d= 2, 
j-r- r> ifrf = 3, 

47r|x — y| 

be the fundamental solution of the Laplace operator. Let then 

S\ ■= J$( ■ ,y)A(y)dr(y), 
V<p := f ^(■ ) yMy)dT(y) l 



Jr 

be the single and double layer potentials, that define solutions of the Laplace equation in 
R d \ T for arbitrary A G #~ 1/2 (T) and ip G H 1/2 (T). Decay at infinity of SX and Vip is 
0(l/r), with the additional assumption that X & H (T) in the two-dimensional case. 
The exterior solution can be represented using Green's Formula (Green's Third Identity): 

u+ = V^u + — Sd u u + = V(^u — Po) — SX, X\=d p u + . (2.4) 

Considering the boundary integral operators 

VA:= 7 ± (5A) K(p:= §( 7 + £V + 7 -£V), (2.5) 

the transmission problem can be equivalently written as the search for (u, A) G if 1 (fi_) x 
H- l ' 2 (r) that satisfy 

(Vu,Vv) n _-<A, 7 v)r =(/,v)n_ + <ft,7f>r Vf G iJ 1 ^-), 

( / i,|7^-/C 7 M)r + ( yU ,VA)r = §A> - £A>)r V^r^fr), 

followed by the integral representation (12.4ft . Note that in the two-dimensional case, 
condition (12. 2p guarantees that A G Hq 1 ^ 2 {T), as can be seen by testing the first equation 
in (12. 6 p with v = 1. An important aspect of the single layer operator V is the fact that it 
is coercive, namely, there exists Cr > such that 

Cr'UW^r < (A,VA> r = ||V(SA)|| 2 d { ™ |^(r) ^ = 3! ^ 

In some forthcoming arguments (but not in the numerical method itself) it will be useful 
to decompose the interior unknown as 

u = 11+ + c c G P (ft_) G := {v G L 2 (fi_) : (v, l) n _ = 0}. (2.8) 

Since /CI = — |, using the decomposition (12. 8 p in the second equation of (I2.6P tested with 
H = 1 yields the formula 

c = J (VA + |( 7 ^ - ft) - AC( 7 u* - A,)) • (2.9) 
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2.2 Discretization with IPDG and BEM 



Consider now a conforming shape regular family of triangulations {Th} of made up 
of triangles/tetrahedra. Let Eh be the set of all edges/faces (in the two/three-dimensional 
cases) of elements of Th and let be the set of interior edges/faces. On each internal 
edge/face we assume that a fixed orientation of the normal vector has been chosen. Using 
this orientation, to each e G we can associate two elements K± G Th such that e = 
K + n K- and that the normal vector on e points from K- to K + . This allows us to 
introduce the following notation for jumps of traces and averages of normal derivatives 

\u] := u\ K+ - u\ K _ id u ii} = \{d v u\ K+ + d y u\ K _). (2.10) 



This notation coincides with the one used in [19] and differs from the one in the unified 
presentation of DG methods in [3] where the jump of a scalar magnitude is a vector 
pointing in the normal direction and the average of vector valued quantities is a scalar 
magnitude. For the purpose of our analysis (and for ease of comparison with [TH]) the 
choice (12.1 Op seems to be the adequate one. Note, however, that the algebraic expressions 
for the IPDG method are the same with the two choices of notation. The triangulation 
Th creates a partition/triangulation Th of the boundary T. 

The discrete spaces that are needed for the method are simply discontinuous piecewise 
polynomial functions: 

V h := {v h :Q^R: v h \ K e F k (K) \/K G T h }, 
A h := {X h : T^R : X h \ e G P fc (e) Ve G T h }. 

For functions that are sufficiently smooth on each element, we can define the IPDG bilinear 
form: 

ODG(wft,«h) : = (V^, V h v h ) Q _ - ^2(id v u h }, lv h ])e 

-Z'EmVhh M)e + £ ^(KJ, M)e. 

eef° ee£° h e 

Here, is the gradient operator applied elementwise, h e is the length/diameter of e G £h 
and a e > is a constant associated to each edge/face. Careful explanation of how to arrive 
at this bilinear form (and why), as well as of the three choices for £ G { — 1,0,1} can be 
found in [3]. A quantity we will need to control is 

cr min := mino- e . 

The coupled IPDG-BEM scheme consists of the non-conforming Galerkin discretization 
of ( 12. 6 p based on the subspaces Vh and A^, and on the substitution of the Dirichlet form 
(Vu, Vu)n_ by aj)c(u,v). The method looks for (uh, Xh) £ V/, x A/, such that 

avG{uh,v h ) - (X h ,v h ) r = (f,v h )n_ + (pi,v h )r Wv h e V h , (2 11) 

(Uh, \u h - Ku h ) T + (ji h , VX h )r = (Hh,iPQ-£Po)r V//^, G A h . 
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All occurrences of angled brackets in (12. lip correspond to L 2 (T) inner products (they are 
not duality products any longer). In comparison with (12. 6p . the trace operator has been 
eliminated, since now Uh and Vh are piecewise smooth functions that can be restricted 
to the boundary, although not with the global trace operator. The additional regularity 
assumed for 0i (for the sake of existence of solution to the original problem, only (3i e 
i/~V 2 (r) i s needed) is justified by its occurrence in the right-hand side of (12. lip . The 
operator /C acts on the restriction of u h to the boundary. Since /C : L 2 (T) — > L 2 (T) is 
bounded (see [H]), this term poses no problem from the practical point of view. 

Algorithmic aspects of this formulation are treated in [19]. Let us just point out a 
detail referred to the two-dimensional case that is not covered in [T9] . If we test the first 

________ 1/2 

equation of (12. lip with % e 1 g Vh, it follows that X h G H (T). This allows for the 
discrete exterior solution 

u+, h :=V(u h -p )-S\ h (2.12) 

to have the right 0(l/r) decaying behavior at infinity. For some future theoretical con- 
siderations, it will be convenient to consider the space 

A ( 0)= f A,n#- 1/2 (r), if d = 2, 

h \ A h , if d = 3. 

Also, parallel to (12. 8p . and only for analytical purposes, it will be convenient to consider 
the decomposition 

u h = u h ^ + c h c h e Po(fi-) %6K:=14nL^4 (2.13) 

Note that 

Ch = ~J (vx h + ±K* - A,) - £K* - A))- (2.14) 

2.3 Main results 

In order to deal in a simpler way with the analysis of the method (12. lip , we introduce a 
global bilinear form: 

B((u h ,\ h ),(v h ,n h )) := a BG (u h ,v h ) - (X h ,v h ) r 

+{\u h - lCu h , Hh)r + (VX h , /i h ) v . 

In principle, this bilinear form is restricted to discrete elements in both components. With 
some abuse of notation, we will allow the exact solution of (12. 6p to be placed in the first 
component of this bilinear form. This can be done assuming some additional regularity for 
this solution and recuperating trace operators for the restrictions of u to V. An important 
property of this method is its consistency: 

B((u h ,X h ),{v h ,Hh)) = B((u,X),{v h ,fi h )) V(u h , nn) e V h x A h . (2.15) 

For analytical purposes we now need to restrict the kind of meshes. We assume that there 
is a neighborhood of V where the triangulation Th is quasi-uniform. This fact will be 
needed for the use of an inverse inequality (see (I3.12p below) on the discrete space A h . 
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This might be restrictive in some practical applications when (3$ and /3\ do not vanish. 
Nevertheless, there are situations when the transmission problem (12.11) (or equivalently 
(12. 6p ) originates from cutting off the computational domain for a problem in free space: in 
this situation the boundary T can be placed apart from the support of the data function 
/ and u is continuous across T (that is (3 = and f3\ = 0); for that case, the local quasi- 
uniformity of the triangulation might not be very restrictive. It has to be emphasized 
though, that at present the inverse inequality plays an important role in the coercivity 
estimates and cannot be easily removed. 

The analysis starts with a quasi-coercivity estimate in terms of a discrete seminorm 
and the bilinear form induced by the operator V (see (12.71) ). For convenience, we write: 

\u h \l:=Y.K l \\iu h \\\l 

The detailed proof of this first result requires the introduction of several elements and the 
proof of some preliminary estimates. This will be dealt with in Section [3j The proof of 
the theorem itself will be given in Section 13.41 

Theorem 2.1 (Coercivity). There exists a constant Cq depending on characteristics of 
the mesh (shape-regularity, local quasi-uniformity near T) and the polynomial degree k 
such that if cr min > a , then 

B((u h , X h ), (u h , X h )) > \\\V h u h f a _ + \\u h \ 2 h + l(X h , VX h ), VK, X h ) eV h x A h 0) . 

Corollary 2.2 (Unique solvability). Under the hypotheses of Theorem \2.1\ the system 
(12. lip has a unique solution. 

Proof. Since (12.111) can be equivalently written as a square system of linear equations, 
only the uniqueness part is needed. If (uh, Xh) G V/, x A/, is a homogeneous solution of 
(12. lip , then testing with Vh = 1 shows that A^ G . Theorem 12.11 and (12. 7p show then 
that Xh = and Uh = c. The decomposition (12.1 3ft - ( 12~T4]) shows finally that Uh = 0. □ 

Convergence estimates will be only given for a solution of the highest regularity. For 
the interior field u, we will consider 

H m {T h ) := J] Hm ( K )> Kt„ := E K*- 
KeT h KeT h 

For the boundary unknown A, we consider the set of edges/faces {Ti, . . . , T M } of the 
polygon/polyhedron T and the broken (but not discrete) Sobolev spaces 

M 

X m (T) :=l[H m (r e ), 

i=i 

endowed with the product norms. 
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Theorem 2.3 (Convergence). Assume that the solution of (12. 6 j) m H k+l {7~h) xX i+1 (r). 
Then, the error of the method of (12. lip can be bounded as 

\\u - u h \\ n _ + ||Vu - V h u h \\n_ + \u h \ h + ||A - Xh\\-i/2,r < Ch k (\u\k+i,r h + |A|fc+i,r) • 
Moreover, for all x £ fi + , 

|it+(x) - u +)fc (x)| < C7/||$(x, ■ )|| 2 , n _ (|«| fc+ i, Th + |A| fc+ i, r ) ■ (2.16) 

As will be made clear by the proof of this result (Section I4.3p . the estimates of Theorem 
12.31 can also be given in more local terms both for u and A, using local meshsizes. That 
part of the analysis is just related to approximation properties and we will not insist on it. 
Note that the bound \uh\h = 0{h k ) takes into account how close Uh gets to be continuous. 

3 Discrete coercivity 
3.1 A geometric construction 

Part of the forthcoming analysis hinges on a class of rigid scaling arguments in order 
to deal with solenoidal vector fields and fractional order Sobolev spaces with different 
scalability properties. We assume the existence of a finite set of reference configurations 
with the following characteristics: 

• In the two-dimensional case e := (0, 1) x {0} and K is a fixed isosceles triangle with 
e as base and such that its two equal angles are less than half the minimum angle 
of the sequence of triangulations Tu- By rotation, translation and scaling, we can 
place an isosceles triangle K e , congruent to K, with base on every e £ and such 
that K e C In order to have unified notation with the three-dimensional case 
we will write /C := {K}. 

• In the three-dimensional case we take a set of pyramids /C := {ifi, . . . , K#} with 
respective bases {ei, . . . , e^}, with C 1R 2 x {0}. We assume that every interior 
face of the triangulation e G S%, with diameter h e , can be extended to be the base 
of a pyramid K e C fi_, such that h~ 1 K e is a rigid motion of one of the elements of 
K. 

Given now any e 6 ^, we consider an invertible affine map F e : M. d — > M. d of the form 

F e ( x ) : = /i e e e x + b e , eje e = i d , 

such that F ( T 1 (i^ e ) G K, and F~ 1 (e) is contained in the base of F~ 1 (i^ e ). Note that in 
the two-dimensional case, this implies that F~ l (e) = e, whereas in the three-dimensional 
case F~ l (e) is a triangle contained in the base of F~ l (K e ). By construction, there exists 
C > (in the two-dimensional case C — 1) such that 

EH* e <Co|MlL WGL 2 (fi_). (3.1) 
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3.2 Fractional order Sobolev norms in reference configurations 

Given a bounded open set OcM d (or the closure of a bounded open set) and s G (0, 1) 
we consider the space H s (0), endowed with the Sobolev-Slobodetskij norm 

IMUo == « + Ko, \v\lo :=[ f fe^ dxdy (3.2) 



oJo i x ~yi 



and the space H s (0), endowed with the norm 



/ ^ 



The negative order spaces are defined by duality, considering the duality representation 
that pivots around L 2 (0), 

H- s {0) := H s (0)', H- s (0) := H s (0)', 

and endowed with the dual norms: 

ii ii (u,v)o .. .. (u,v)o . , 

IMI-s,o := SU P Ti — f] ) IMI-s,o,~ := SU P fi — f] • (3.4) 

For open polyhedral surfaces (polygonal curves), we can consider the spaces H S (S) and 
H S (S) for — 1 < s < 1. For closed polyhedral surfaces T, the spaces H S (T) for < s < 1 
are defined by parametrization and H~ S (T) := H S (T)', pivotal to £ 2 (r). 

Since for e E (0,1/2) the trace operator H 1 ~ £ (K) — > H 1 ^ 2 ' £ (dK) is bounded and 
surjective, it has a bounded right inverse L £ and we can bound, with a constant C® > 
that depends only on e, 

\\V(LM\- £ ,k<C®M 1/2 _ £:9R VtpeH^idK), Ket (3.5) 

Since H S (K) = H S (K) for < s < 1/2, with equivalent but different norms, for e G 
(—1/2, 1/2), we can identify constants C® > depending only on e such that 



,K,~^ C ?\M-e,K VveH-%K), Ket (3.6) 



The third group of inequalities requires some additional work in the three-dimensional 

case. 

Lemma 3.1. For e G (0, 1/2), there exists Cf> > such that 

Mi/2-e, 9 K < C?|M|i/2- £ ,a~ V</? G H 1 ^ 2 ~ £ (D), K G Z, (3.7) 

where D = e if d = 2 and D is any triangle contained in the base of K if d = 3, and <p is 
the extension by zero of <p. 
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Proof. The first part of the proof is only required for the three-dimensional case. Let e 
be the base of K and D be any triangle contained in e. Then 

1 2 I 1 2 , o f I f ^ X I / \| 2 



< + 2 / (2W * ) My) |My 

J D \ Jdist(y,dD) ' 

|2 , 2vr / |v,(y)| 2 



and therefore 



£ l < D + T/ D d£t(£aD)5 dy ' 



where the key fact is that the constant does not depend on D. In the second step, 
we apply that extension-by-zero is a bounded operator from H s (e) to H s (dK) and take 
s = l/2-e. □ 

Lemma 3.2. For e G (0, 1/2), t/iere exists C® > st/c/i that 

||v • n||_ 1/2+£iD < C®||v|| £) £ Vv G H £ (K) d such that divv = (3.8) 

where D = e if d = 2 and Z) any triangle contained in the base of K if d = 3. 

Proof. Let 9? G H 1 ^ 2 ~ £ (D), lp be the extension by zero of ip to the rest of dK, and 
u := L e <£>, where L e is the trace lifting of (13.51) . Since the hypotheses on v imply that 
v G H(div, K), we can apply Green's Theorem and prove that 



(v-n,<p) D = (v • n,ip) dR = (divv,u) R + (v, V 



u 



K 



< \WI,kW Vu \Ur,~ ^ c?IMI e ,*llv«IL i * 

< ^®C?l|v||^||^|| 1/2 _ Ei aK < C? Cf Cf ||v||^||^|| 1/2 _ £j ^, 

where we have applied (13.61) , (13.61) and Lemma 13.11 The proof is now a direct consequence 
of the definition of the H~ l / 2+e (D) norm. □ 

3.3 Two key lemmas 

We first divide the set of edges as 

C iP := {e € El : dist(e, T) < C} and := £° h \ (3.9) 

If C above is large enough (or the triangulation is refined enough), we can fit 

U {K e : e G C O int C (3.10) 

Since we are assuming that the triangulation is shape-regular and quasi-uniform near the 
boundary T, then the partition is quasi-uniform and the diameter of the elements of 
is equivalent to 

/istrip := max{/i e : e G ££ tnp U T h }. (3-H) 
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Therefore, we have an inverse inequality for elements of the discrete space A^: 

K tlip \\X h \\- 1/2+e , T < C™\\\ h \\„ 1/2tT V\ h EA h , e G [0, 1/2]. (3.12) 

An additional ingredient for the proofs below is related to continuity of layer potentials, 
namely, we can choose C® for every e G [0, 1/2] such that 

||V(SA)|| £Ant < C?||A||_ 1/2 , r VA G H^\Y) (3.13) 

and 

||V(<SA)|U_ < C©||A||_ 1/2+£ , r VA G H-^+%T). (3.14) 

Note that the interior regularity bound (I3.13P can be proved directly by bounding the 
fundamental solution, since the distance between Q- mt and T is positive. On the other 
hand, (I3.14p is a well-known regularity result of layer potentials (see [TT]). 

Lemma 3.3. For all e G (0, 1/2), there exists C® > such that 

J2 h d - 2 \\d u (u* o F e )|| 2 1/2+£>F - 1(e) < C®\\Vu% d Vu* = SX h , X h G Af . 

Proof. Note that Am* = in K e for all e, since K e C fi_. Since F e is a dilation and a 
rigid transformation, it is clear that A(u* o F e ) = in K. Therefore, by Lemma [3.21 we 
can bound 

R(u* o F e )||_ 1/2+£jF - 1(e) < C®\\V(u* o F e )|| a Ve G £° h . (3.15) 

A simple change of variables proves that 

\\V{u*oF e )\\l k = \\V(u*oF e )\\l + \V(u*oF e )\l R 

= e d (l|Vull^ + /aV^ e )- (3-16) 

Before adding over all edges, notice that the estimate (13. ID and the definitions of the 
subsets of edges (13.91) and property (13. 10[) imply that 

\ V \e,K e < C o\v\e,Q iat and \ V \e,K e ^ C oMe,U.- 

e£E int eee s h tIip 
Hence, (13.161) . the estimates (13. 131) - (13. 141) and the inverse inequality (13 . 1 2[) imply that 
h d e - 2 \\V(u* o F e )\\l R < C (\\Vu*\\l_ + h^\Vu% nint + h% ip \Vu% n _) 

< a(||Vn1|L + (C'®) 2 (/ i 2£ + Cf- 2£ )||A,|| 2 _ 1/2)r ). (3.17) 
The result is then a straightforward consequence of ( 13 . 1 5j) . (13. 1T|) and (12. 7\i . □ 
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Lemma 3.4. For all e G (0, 1/2), there exists C® > such that 

£ ht 2 \\lu h o F e ]\\l /2 ^ F - He) ^ < C®\u h \ 2 h \fu h G V h . 

Proof. Some preparations are first needed in the three-dimensional case. Let D be a 
fixed triangle in 1R 2 x {0} and let D e := F~ 1 {e). The affine maps G e : R 2 -> R 2 that 
transform D to D e are uniformly bounded and have uniformly bounded inverses, because 
the triangles D e are shape regular and their area is of order one. Therefore, we can bound 

\\<p\\.,3 e ,~^ B A\<P oG e\\.£,~ VcpeH s {D e ), 0<s<l. (3.18) 

Using ( I3.18P with s = 1/2 — e, a finite dimension argument on D, and the change of 
variables G e , we can prove that there exists C® such that 

INIW^ < C ?M% V« G P fc . (3.19) 

Note that this result is straightforward in the two-dimensional case since then for all e, 
D e = e. Finally 

h d ~ 2 1| [u h o F e j |i; /2 _ £jF - 1(e)i . < ht 2 C®\\ [u h o FJ ||| e = CfK 1 1| M || 2 . 

□ 

3.4 Proof of Theorem I2J1 

A standard argument in DG analysis (see [3] for example) can be used to prove that 

ODG(ufc,Wh) > l|V A Mfc||n_ - |1 + ^l(c*^l|V/ l M/ l || 2 7 _ + {A8)~ l \u h \l) + o- min \u h \l, (3.20) 
where C* is a positive constant that allows us to bound 

MKIIe < C*IMlL Vv h eV h . (3.21) 

K€T h £ h €eCdK 

Taking 5 = 1/(4C*|1 + £|) in (13.201) for £ ^ —1, we can present all cases together with the 
estimate 

QdgK,0 > |I|V^|| 2 _ + (a min - |l-K| 2 a)KI' Vu h eV h . (3.22) 

We now turn our attention to the full bilinear form. Let G (satisfying the additional 
constraint j r X h = in the two-dimensional case). Let then u* := SXh- By well-known 
results on boundary integral operators it follows that 

±\ h - K'Xu = -dU\ \ h = d;u* - d+u*, (VX h , X h )r = || V«*|| Kd (3.23) 
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(recall ( 12. 7p ). Therefore 

B((u h , X h ), (u h , X h )) = a DG (u h , u h ) - (X h , u h ) T + (u h , \X h - /C%) r + (VX h , X h ) r 

= a DG (u h ,u h ) - (d~u* - d+u*,u h ) r - (d+u*,u h )r + ||Vw*||jjj, d 

= a DG (u h ,u h ) - (VV, VhU h )n_ + IIWIIr* + ^ ({{<9„w*}}, [u h j) e , 

where in the last step we have applied Green's formula element-wise, taking advantage 
of the fact that Au* = 0. Using (I3.22f) and the fact that d u u* is single-valued across 
inter-element faces, it follows that 

B((u h , X h ), (u h , \ h )) > J||V^||L+|||Vw*||^ 

+ Kin " |1 + £| 2 O0Klh. + ^{dyU*, M) e . (3.24) 

Finally, since F e are rigid motions composed with dilations, we can easily apply Lemmas 
13.31 and 13.41 to bound 

E<^*,M) e = ^^ 2 (9,KoF e ),[« h oF e ]) F - 1(c) 
eee° h e ee° h 

< E hd e~ 2 W d «( U * ° ^)||_x/ 2+e> ^(e)IIK O F e j\\ 1/2 - £ , F -^ 

n® 

< c®5\\w\\l d + ^\u h \l 

Taking now 5 = C®/4 and inserting this bound in (13.241) . it follows that 

B((u h , X h ), K, X h )) > l\\V h u h \\l_ + \\\Vu*\\ 2 Rd + (a min - |1 + e| 2 a - C®/Cf)\u h \ 2 h . 
The result now follows by simply taking <r m i n > ~ + |1 + ^| 2 C* + C®/C® and using (12.71) . 

4 Convergence analysis 
4.1 Quasi-optimality 

We consider the following discrete seminorm, defined in H x {%) x H-^ 2 (T) 

\i(u,X)t\l;= ||V,n||L + |n|^+||A|| 2 _ lAr , 
as well as the stronger norm 



ik«. mist ■■= IK mi + E k'hi + E + E h 

defined in H 2 {T h ) x L 2 (T). 



2 L.|2 

IL 

eer h e gr h KeT h 



K\ U \%Ki 
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Proposition 4.1 (Continuity of the bilinear form). There exists G\ > such that 

\B((u,\),(v h ,fi h ))\ < CMu^ltv&fariU l^l^f^xlf^ 
Proof. Using a local trace inequality 

£ h M\l < h K \Hl K < (\\v\\l + hl\\Vv\\%) VveH\K), (4.1) 

£ h 3eCdK 

and applying ( 13.211) to the components of Vh v hi it is simple to obtain the bound 

1/2 



v h )\ < \\V h u\\ Q _\\V h v h U-+CH 2 \\V h u\\l_ +Y, h 

V KeT h 



kHIk I \Vh\h 



+Cl /2 \u\ h \\V h v h \\ n _ + a max \u\ h \v h \ h . (4.2) 
Also, using a discrete Poincare-Friedrichs inequality (see jl] and [2]), we can bound 

Hlln_ < CW(||V^||L + \v h \ 2 h ) V2 Vv h e V£, (4.3) 
which, together with (I3.2ip allows us to bound 

1/2 f \ 1/2 

2 



l(A,« fc >r| < fE^IIAIlA fe^K 

\eer h / \eer h 

< c^CvAy^K^AZ] (llv^HL + Kg) 

\eer h J 



h 
1/2 



1/2 

(4.4) 



The third term in the bilinear form is bounded using the inverse inequality (13.121) (recall 
(13. lip for the definition of the local meshsize), proving that 

\{\u - Ku,n h ) T \ < \\\X- /C||z,2(rHL 2 (r)||«||r||A*/i||r 

< Ci7 2 ||iX-/C|| L2(rHL2(r) (y,K X \\u\\1\ IMI-W- (4.5) 

Ver h J 



The result is now a direct consequence of (14.21) . (I4.4p . ( I4.5p . the continuity of V : 
F _1 / 2 (r) ->• H 1/2 (T) and the definition of the broken norms. □ 

Let Il/j : L 2 {VL_) — > Vh and Ph : £ 2 (r) — > A h be the orthogonal projections onto Vh 
and A h respectively. Since constant elements are in Vh, we can decompose (using (12. 8p 
and (gJSD ) 

u h - Yl h u = u h> ± - + c h -c. (4.6) 

Also, using the fact that constant functions are in A^, we can show (see the first equations 
in (EED and (124H ) that 

f P h \ = I \= I f + f A = / X h . (4.7) 
Jr Jr Jn^ Jr Jr 

(All of them vanish in the two-dimensional case because of (12. 2p .) 
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Proposition 4.2 (Quasi-optimality). Let (u,X) and K, A/J be the respective solutions 
of ( 12"H]) and (12. lip and assume that (u,X) G H 2 (Th) X L 2 (T). T/ien i/iere exists C 2 > 

IK - w, A fc - A) ||U < C 2 \l(u - U h u, X - P h X)\\ hjSt . 

Proof. Note that by KT} . X h - P fc A G A^ 0) . We can then apply Theorem O and to 
bound 

IK - u hU , x h - p h x)f h = IK, - n h u*, x h - p h x)f h 

< 4max{l,C r } J B(K iJc - 11^, X h - P h X), K,* - IL^, X h - P h X)) 

= 4max{l, C T }B{(u - U h u, X h - P h X), K* - U h u*, X h - P h X)), (4.8) 

by the consistency of the method (12.151) . the decomposition (14. 6 p and the fact that 

B((c - c h , 0), (u ft) * - U h u+, X h - P h X)) = (c - c A ) J{X h - P h X) = 

(see (14 .7p and recall that /CI = — ~). Since w^* — H^w* G V^, we can apply Proposition 
14.11 and the result is a direct consequence of (14. 8 p . □ 



4.2 Estimate of the interior average 

Because we are dealing with the pure transmission problem (no boundary conditions), 
the discrete seminorm does not take into account the average of the component u of the 
solution in a. This is very clear from the coercivity estimate (Theorem 12. ip . which is 
written in terms of that seminorm and, as such, is not able to estimate the error 



Ch 



a 



(u - u h ) 



(see decompositions (12. 8p and (I2.13P ). We start the analysis of this term with a lemma 
that combines some of the arguments of the previous sections. 

Lemma 4.3. Let £ := | — K}\. Then there exists C% > such that 

\(Z,v h ) r \ < C 3 (||V^||L + \v h \l) V2 \fv h G V*. (4.9) 
Proof. Note that £ = 1 — d~Sl. We can then bound 

\(d-Sl,v h ) r \ < C A (\\V h v h \\L + \vh\ 2 h ) 1/2 Vv h G V h , (4.10) 



by proceeding as in the proof of Theorem 12.11 (see Section |3~4|) . Note that when we aim 
to apply Lemma 13.31 in the two-dimensional case X h = 1 £ A h 0) . However, we can still 
bound as in (13.171) 

£ ^- 2 ||V((51) o F e )\\l R < C (||V(S1)||L + h^\V(Sl)\l n 
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which is enough for our current needs. In order to estimate the remaining term, we 
consider the solution of the interior Neumann problem 



—Am = c in Q , d v u = 1, with c := 



iri 



and note that classical Sobolev regularity results can be invoked to show that u G 
H 3 / 2+£ (Q_) for some e > 0. We then apply integration by parts element-by-element 
to obtain 

(1, v h ) T = -J2 ( d »u, Me + (W, V h v h ) n _ - c(l, v h ) Q _ . (4.11) 



The first term is bounded using 



~ y^(di/u, jvhDe < ( ^2 h e \\d„u\\f) \v h \ h < (^^2 h K \\Vu\\ 2 dI< ) \v h \ h 



K€T h 

1/2 : 



< Cf(Y,\\Vu\\ 2 K + h%\Vu\l K ) \v h \ h , (4.12) 



K€T h 

where the last inequality follows from a scaling argument and the trace inequality in a 
fixed reference element. Since Vh G Ll(Q_), the result follows from the inequalities (I4.10p . 
fr4~TT|) and KT2\> . □ 

Proposition 4.4. There exists C5 > such that 

\c-c h \< C\\{u-U h u,X-P h X)\l h>st 



Proof. Using the integral formulas for c and Ch (that is, ( 12 .9p and (I2.14p ). we can easily 
bound 

|r| \c-c h \ < || V(A — A h )||i /a>r +||e||r||7W*-n /l M*||r + K^«* J /,-n h u*)r| 

< l|V||£r-va(rHfli/a(r)P ~ K\\-x/2,r + ||£||r||7« - n Hlr + - n ft it*) r | 

where £ := | — /C*l. Note that it^* — TlhU+ G and we can thus apply Lemma [O] to 
bound 

|<£, - n ft «*) r | < QIIK,/, - n h it*, A ft - p A A)[ h = IK - n A «, x h - P h x)\\ h . 

The result is now a consequence of the quasi-optimality bound (Proposition 14.2ft and the 
definition of the discrete norms. □ 



4.3 Proof of Theorem 12^1 

By simple approximation properties (using scaling arguments and the Bramble-Hilbert 
lemma), we can bound 

\\(u-U h u,X-P h X)\\ h>st <C 6 ( Yl h 2 Ml + i,K+Y, h e k+1 \Ml+i,e) 1/2 (4-13) 

KeT h eer h 
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This inequality and Proposition 14.21 provide the bound 

1 /2 

||V« - V h « fc ||n_ + \u h \ h + ||A - A fc ||_ 1/2 , r < C 7 ( ^ ^MjUi,* + S ^Wfc+i.e) ■ 

For the bound of the L 2 (f2_) error we use the decompositions (12.81) and ( 12.1 3ft and the 
discrete Poincare-Friedrichs inequality (14. 3 p to prove that 

\\U h u - u h \\n_ < \c - c h \ |fi_| 1/2 + HriftW* - u h ,*\\n_ 

< \ c -c h \ \n^\ 1/2 + c PF i\(u h u,-u h ^o)\i h 

Propositions 14.21 and I4.4[ and the approximation estimate (14.131) then provide the L 2 (fi_)- 
estimate. Note that this estimate is suboptimal (some of the IP methods included in this 
analysis have no super convergence properties), but that it is needed as a complement of 
previously given error bounds, since ||| • is only a seminorm. 

Subtracting the representation formula (12 .4p from the discrete approximation for the 
exterior solution fl 2 . 1 2 j) . and using the definition of the layer potentials given in Section 
12.11 we can bound 

|u +ifc (x)-ii+(x)| < ||$(x, • )||i/ 2 ,r||A- Aj_ 1/2) r + ||$(x, • )|| r ||7« - n Hlr 

+ |(^$(x, • ),U h u-u h ) r \. (4.14) 

For fixed x G fi+, $(x, • ) G C°°(f2_) is harmonic. Using element-by-element integration 
by parts and a scaling argument (14. lft . we can write 

|(^$(x, -),v h ) r \ = |-£<$,$(x, .),N)e + (V$(x, O.Vfc^n. 

1 /2 

< ( E ^-||V$(x, • )||^) KL+ ||V$(x, ■ )||nJ|V^|| n _ 

< max{CV 2 , h 2 CU 2 + l}||$(x, • )|| 2>n _|K0)|U V« fc G V fc . 

The bound ( 12. 16ft is then a consequence of this latter inequality applied to Vh = TlhU — Uh, 
( I4.14p and the above estimates for the error. 
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